RR10_test_cerebelum_ST_clustering

0.1 Loading packages, creating seurat object

0.1.1 Set paths

Code

# Set parent dir for visium experiments containing folders of spaceranger output with names
visium.dir <- "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RRRM2_RR10_test/RR-10/Cerebelum"

# Path for where the plots will be exported to
export_path = "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RRRM2_RR10_test/RR-10/out_cereb"

0.1.2 Create infoTable and create Seurat object

Code

samples <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'filtered_feature_bc_matrix.h5')
spotfiles <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'tissue_positions.csv')
imgs <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'tissue_hires_image.png')
json <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'scalefactors_json.json')

Code

section.name <- samples
section.name <- gsub(paste0(visium.dir, "/"),"", gsub("/filtered_feature_bc_matrix.h5","",section.name))

Code

infoTable <- data.frame(section.name, samples, spotfiles, imgs, json, stringsAsFactors = FALSE)

1 Data Quality Control

1.1 Number of Genes per Spot

Code

VlnPlot(se_brain, features = "nFeature_RNA", ncol = 1, group.by = "section.name", pt.size = 0)

Code

ST.FeaturePlot(se_brain, features = "nFeature_RNA", dark.theme = F, cols = c("lightyellow", "red", "dark red"), pt.size = 1.2, ncol = 4)

1.2 Number of UMIs per Spot

Code

VlnPlot(se_brain, features = "nCount_RNA", ncol = 1, group.by = "section.name", pt.size = 0)

Code

ST.FeaturePlot(se_brain, features = "nCount_RNA", dark.theme = F, cols = c("lightyellow", "red", "dark red"), pt.size = 1.2, ncol = 4)

1.2.1 VlnPlot - UMI/Genes/mtGenes/RiboGenes

Code

se_brain$percent.mt <- PercentageFeatureSet(se_brain, pattern = "^mt")
se_brain$percent.ribo <- PercentageFeatureSet(se_brain, pattern = "^Rpl|^Rps")
VlnPlot(se_brain, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.ribo"), ncol = 2, pt.size = 0)

Run PCA here on the non-filtered data and then plot it

1.3 Ploting of the non-filtered data

Code

p1 <- DimPlot(se_brain_merged_notfilt, group.by = "section.name")+ ggtitle("Plot of the non-filtered data")
p2 <- DimPlot(se_brain_merged_notfilt)

p1|p2

Code

coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A")

DimPlot(se_brain_merged_notfilt, cols = coldef)

1.4 Filtering step

Code

enids <- read.table(file = "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RR3_test/rr3_ref/mm10_genes.tsv", header = T, stringsAsFactors = T)
enids <- data.frame(apply(enids, 2, as.character), stringsAsFactors = F)
rownames(enids) <- enids$gene_id

Code

se_brain <- SubsetSTData(se_brain, expression = nFeature_RNA > 200)

1.4.1 mitochondrial, ribosomal genes and Malat1 filtered out

1.4.2 export Staffli object

Code

staffli <- se_filtered@tools$Staffli

Code

se_brain_merged <- merge(se_section_list[[1]], y = c(se_section_list[2:length(se_section_list)]), merge.data = T)
VariableFeatures(se_brain_merged) <- features

Plot data before and after harmony

1.5 PCA + Harmony

1.5.1 PCA

Code

se_brain_merged <- RunPCA(se_brain_merged)
se_brain_merged <- RunUMAP(se_brain_merged, dims = 1:35, reduction = "pca")

Code

p3 <- DimPlot(se_brain_merged, group.by = "section.name")+ ggtitle("Plot of the filtered data before integration with Harmony")
p4 <- DimPlot(se_brain_merged)

p3|p4

Code

coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A")

DimPlot(se_brain_merged, cols = coldef)

1.5.2 Harmony

1.6 Clustering

Code

p7 <- DimPlot(se_merged_harmony, group.by = "section.name")+ ggtitle("Plot of the filtered data after integration with Harmony")
p8 <- DimPlot(se_merged_harmony)

p7|p8

Code

coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A")

DimPlot(se_merged_harmony, cols = coldef)

Code

se_merged_harmony@tools$Staffli <- staffli

Code

ST.FeaturePlot(se_merged_harmony, features = "seurat_clusters", dark.theme = F, pt.size = 2.3, show.sb = F, ncol = 4, cols = coldef)

Code

#top5_2 <- markers_RR10_2 %>%
top5_2 <- markers_RR10 %>%
  group_by(cluster) %>%
  filter(p_val_adj < 0.01) %>%
  filter(row_number() %in% 1:5)

d2 <- DotPlot(se_merged_harmony, features = unique(top5_2$gene) %>% rev()) + 
  coord_flip() + 
  scale_colour_gradientn(colours = RColorBrewer::brewer.pal(n = 11, name = "RdBu") %>% rev()) + labs(y = "cluster")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

Code

d2